{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9e33c919-1309-41a1-807b-f0406574dbb5",
   "metadata": {},
   "source": [
    "# Ligand Parameterization\n",
    "## Objective\n",
    "\n",
    "The objective of this tutorial is to do a ligand parameterization. Molecular mechanics forcefields need to be told how to treat each atom via a set of parameters. If there is a molecule (residue) in your system that your forcefield of choice does not already have parameters, you will need to build these yourself. In this example, we will see how QMzyme can automate the calculations required to parameterize a ligand in line with the RESP procedure (J. Phys. Chem. 1993, 97, 40, 10269–1028). \n",
    "\n",
    "This workflow allows you to:\n",
    "\n",
    "- Prepare ligand for ligand parameterization.\n",
    "\n",
    "In this specific example, we are using ketosteroid isomerase (KSI) as the model system. The structure for KSI is obtained from the PDB [1OH0](https://doi.org/10.2210/pdb1OH0/pdb) and MM-minimized prior to this tutorial. We will utilize it to create an input file for ORCA for geometry optimization and frequency calculation.\n",
    "\n",
    "## Classes used in this example\n",
    "\n",
    "- [GenerateModel](https://qmzyme.readthedocs.io/en/latest/API/QMzyme.GenerateModel.html)\n",
    "- [QM_Method](https://qmzyme.readthedocs.io/en/latest/API/QMzyme.CalculateModel.html#qm-treatment)\n",
    "\n",
    "## Required Files\n",
    "\n",
    "To start, you will need:\n",
    "\n",
    "- A fully prepped and protonated PDB.\n",
    "  \n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c90413f-1995-47eb-bd4a-f0b4f6efb2c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here are the necesary imports for this tutorial!\n",
    "\n",
    "import QMzyme\n",
    "from QMzyme.data import PDB"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "428ce54d-3d0c-43f9-91f6-a7565ca92c20",
   "metadata": {},
   "source": [
    "## Initialize Model\n",
    "\n",
    "To begin our example workflow, we need to start by initializing QMzymeModel using `GenerateModel()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "c3f2b5c3-a90d-4291-8393-c92a738d786f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Charge information not present. QMzyme will try to guess region charges based on residue names consistent with AMBER naming conventions (i.e., aspartate: ASP --> Charge: -1, aspartic acid: ASH --> Charge: 0.). See QMzyme.data.residue_charges for the full set.\n",
      "\n",
      "\tNonconventional Residues Found\n",
      "\t------------------------------\n",
      "\tEQU --> Charge: UNK, defaulting to 0\n",
      "\n",
      "You can update charge information for nonconventional residues by running \n",
      "\t>>>QMzyme.data.residue_charges.update({'3LETTER_RESNAME':INTEGER_CHARGE}). \n",
      "Note your changes will not be stored after you exit your session. It is recommended to only alter the residue_charges dictionary. If you alter the protein_residues dictionary instead that could cause unintended bugs in other modules (TruncationSchemes).\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# We first initilize QMzymeModel by loading structure file.\n",
    "# The initial pdb file should be prepared prior (hydrogens must be present).\n",
    "model = QMzyme.GenerateModel(PDB)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10ee63fd-d748-40b4-b7af-c37e0cb9d0f7",
   "metadata": {},
   "source": [
    "## Add Ligand Charge\n",
    "\n",
    "When loading the PDB file, we quickly run into an issue: a nonconventional residue is found. This simply means that there is a non-AMBER protein residue with an unknown charge. This can be resolved using `residue_charges.update()`, where we can manually update the charge of an unknown residue. In this case, the charge of equilenin (EQU) is -1 and we will update it using `residue_charges.update()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "55b1bfd3-b8e0-4e88-9672-81ae8096ca7f",
   "metadata": {},
   "outputs": [],
   "source": [
    "QMzyme.data.residue_charges.update({'EQU': -1})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62fe2630-8b6f-41fc-b01c-ab3721c097b2",
   "metadata": {},
   "source": [
    "## Designate Ligand Region \n",
    "\n",
    "In QMzyme, `set_region()` is used to create QMzymeRegion, a Python class object that is the collection of QMzymeAtom objects. To select protein amino acids directly, we use MDAnalysis nomenclature for selecting our residues of interest!<sup>1,2</sup>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "440699e9-b9c2-419a-be91-21e7f35f80cf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ligand region:  <QMzymeRegion EQU contains 37 atom(s) and 1 residue(s)>\n"
     ]
    }
   ],
   "source": [
    "# We'll select the region of interest.\n",
    "model.set_region(name='EQU', selection='resid 263 and resname EQU')\n",
    "print(\"Ligand region: \", model.EQU)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "579483c9-7d49-4888-a755-375c5bde88bc",
   "metadata": {},
   "source": [
    "## Build the QM Method and Assign to Region\n",
    "\n",
    "Ultimately, QMzymeModel is a collection of QMzymeRegon and a QM method for calculation. This can be achieved by creating a QM method and assigning QM method to the region. You can find more about different QM methods and their usages at:\n",
    "\n",
    "- [QM_Method](https://qmzyme.readthedocs.io/en/latest/Examples/QM-only%20Calculation.html)\n",
    "\n",
    "For the purpose of this example, we will forgo geometry optimization and only perform a single point energy. Calculation and population analysis at the level used in the original RESP procedure (J. Phys. Chem. 1993, 97, 40, 10269–1028).\n",
    "\n",
    "After creating QM method of choice, we can assign the QM method to a QMzymeRegion using `assign_to_region()`. This will make a complete set of QMzymeModel, which can be used to now write a QM input file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "392ab132-4158-4e76-9896-70b56201b76c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "QMzymeRegion EQU has an estimated charge of -1.\n"
     ]
    }
   ],
   "source": [
    "# We first create two QM method, each with desired level of theory.\n",
    "qm_method = QMzyme.QM_Method(\n",
    "    basis_set='6-31G*', \n",
    "    functional='HF', \n",
    "    qm_input='SCF=Tight Pop=MK IOp(6/33=2,6/42=6,6/43=20)', \n",
    "    program='gaussian'\n",
    ")\n",
    "\n",
    "# Now we assign this method to our QMzymeRegion EQU\n",
    "qm_method.assign_to_region(region=model.EQU)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3b093f7-f604-4376-8ff2-ad31f22f3522",
   "metadata": {},
   "source": [
    "## Write Calculation Input File\n",
    "\n",
    "The final step now is to create the calculation input file. This can be done by using `model.write_input()` method, which will create a directory named QCALC and save the input file there. When this method is called, the `store_pickle()` method is triggered, and you will have the corresponding .pkl file for your QMzymeModel. This is especially useful because after your calculation is completed, you can then load the calculation results back into the corresponding QMzymeRegion using the QMzymeRegion method `store_calculation_results()`, passing it the calculation file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "933b4fb0-333d-43c9-b280-e1a21d8935a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# QMzyme will know we only have one region with a calculation method set, \n",
    "# so it will logically create the input file for that scenario\n",
    "model.write_input('EQU_resp')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed2feab9-b61b-4716-bb6d-923575a3d303",
   "metadata": {},
   "source": [
    "## References Cited\n",
    "1) N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein. MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. J. Comput. Chem. 32 (2011), 2319–2327. doi:10.1002/jcc.21787\n",
    "2) R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler, D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein. MDAnalysis: A Python package for the rapid analysis of molecular dynamics simulations. In S. Benthall and S. Rostrup, editors, Proceedings of the 15th Python in Science Conference, pages 98-105, Austin, TX, 2016. SciPy. doi:10.25080/Majora-629e541a-00e"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}